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Abstract 

The aim of this paper is a short survey of models and methods that de- 
veloped by the authors. These models and methods are used to optimize 
general networks with nonlinear non-convex restrictions and objectives 
possessing mixed continuous-discrete optimization variables. There are 
discussed the problem formulations and solution methods for simulation, 
optimization, sensitivity and stability analysis for flow with nonlinear po- 
tential in general networks. These problems and the developed methods 
and programs have industrial application e.g. by gas networks. 



Contents 

1 Introduction Q 

2 Mathematical model and problem formulation la 

2.1 Passive and active elements [2| 

2.2 Steady state continuous-discrete optimization problem in general 
network [3] 



Methods of the discrete-continuous nonlinear network optimiza- 
tion 

3.1 Simple examples 

3.2 The continuous dominant problems 

3.3 A method for solving the continuous problems in general networks 

3.4 Modification of the branch-and-bound method 

3.5 Another procedure for searching an optimal integer-feasible solution 



4 Risk, sensitivity and stability analysis for reliability study @ 

5 Conclusions 

*Dr. L. A. Ostromuhov, Wingas, Friedrich-Ebert-Str 160, D-34119 Kassel, Germany, 
t e-mail: leonid. ostromuhov&wingas. de 

* Presented on 16th IMACS World Congress 2000 on Scientific Computation, Applied Math- 
ematics and Simulation in Lausanne, Switzerland. — IMACS World Congresses on Computa- 
tional and Applied Mathematics and Applications in Science and Engineering 



1 



1 Introduction 



There are 2 approaches to describe flow in real and virtual networks presently. 
The first approach is based on a physical description of flow laws. Conservation 
laws of mass, impulse and energy resulting in Ohm law for electrical current or 
Bernoulli law for fluid flow, and thermic state equation are used there. This 
approach led to physical and technical modelling and simulation of real existing 
systems. Dynamic optimization problem is not solved. Instead of dynamic 
optimization, 'try and error' method based on simulation was proposed. 

The second approach is developed more for economics. It is based on graph 
theory and special methods of linear programming. The flow conservation law 
represents a main part of restrictions in a mathematical model there. Such 
mathematical problems as maximal flow in network and minimal cost network 
flow could be mentioned there. In their dual problems, linear potential is ap- 
peared. Linear potential is connected with flow by means of linear analogue of 
Ohm law for passive electric chains. 

However, there exist technical systems within which the connection between 
flow and potential is not linear. The objective function there depends not only 
either on flow or on potential, but on both of them. There are non-convex 
and non-smooth functions among constraints and objectives. The pipe net- 
works, energy systems and electrical chains with nonlinear associated facilities 
are examples for such large - scale technical systems. There are some working 
approaches to optimize such systems [5], [11] . [8] but the optimization theory 
based on dynamic models is not developed for general networks with nonlinear 
connection between flow and potential. 

The aim of this paper is a short survey of models and of developed by 
me methods using to optimize general networks with nonlinear non-convex 
restrictions and objectives possessing mixed continuous-discrete optimization 
variables. There are discussed the problem formulations and solution methods 
for simulation, large-scale nonlinear discrete-continuous network optimization, 
sensitivity and stability analysis and reliability study for flow with nonlinear 
potential in general networks. The developed methods belong to mixed integer 
nonnlinear programming (MINLP) in networks. They use graph theory, inte- 
ger and nonnlinear programming as their bases corresponding to the nature of 
problems. 

The mentioned above problems find an industrial application. The devel- 
oped methods are realized in software called ACCORD which is in operation in 
different countries. 



2 Mathematical model and problem formulation 

2.1 Passive and active elements 

A passive element is a two-pole that does not need an additional energy supply 
to work. As simple examples, the resistors and pipe could be mention. A 
dynamic model of a pipe could be written as a system of algebraic and partial 
differential equations (PDE) for state variable vector (p,T, p, v, h) T : 

P 

(1) p = p (p, h) := , where z — z (p, T) 
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where c p - specific heat by constant pressure [J /kg K] ; M - mass flowrate 
[kg/s]; g - gravitational acceleration [to/s 2 ]; h - specific ehthalpy [J/kg]; F - 
objective function; H - altitude [to]; p - pressure [iV/m 2 ]; R - individual gas 
constant [J /kg K]; p - specific mass (density) [kg/m 3 ]; S - cross section [to 2 ]; 
T - temperature [K]; t - time coordinate [s]; v - flow velocity [m/s]; x - length 
coordinate [to]; z - compressibility factor [— ]; 51 - heat flow per length [W/m]. 
The equations express respectively thermic and calorific state equations and 
conservation laws of mass, impulse, and energy. The PDE is of hyperbolic type. 
For a practical simulation, a mixing from initial and boundary conditions is 
used. 

For a lot of considerations, a steady state model of a passive element is 
precise enough. The model is algebraic and could be obtain after simplification 
from as follows: 



(4) 



q(x) = Const, i < x < k 



(5) 



f(pi,Pk, q) = 0, where Pi = p \ x =i, Pk = P \x=k 



As example of active elements, an amplifier, a pump, and a compressor could 
be mentioned. Their common feature is that every of which have a plane simply 
connected operational range [9]. This range must be taken as a restriction in 
the model. Unfortunately, the range forms a non-convex domain mostly. 

In some cases, in particular for hydraulic networks, active elements react 
much quicker than passive elements. Then steady state models could be used for 
active elements by network modelling. These steady state models are algebraic. 

Together with flow balance equations, the equations and inequalities describ- 
ing models of active and passive elements form a part of restrictions by network 
modelling. 



2.2 Steady state continuous-discrete optimization prob- 
lem in general network 

Let be given a network G — (V, E) with a set of nodes V and a set of edges E. 
Let qik is a flow in the edge (i, k) and Qi is supply or demand called intensity 
of node i. It means that 



(0) 



k£V 



qik +Qi=0, ie V, 



(7) 



qik 



-q k i, (i,k) e E. 
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Suppose that there are given families of functional dependencies on flow and 
potentials: 



(8) fd ik (Pi,Pk,qik,Cik) = 0, (i,k)eE; 

(9) dik € {1, Nik} > (i,k)eE. 

Here Cjfc is a vector of continuous parameters (coefficients), and dik is a discrete 
parameter of the edge (i, k). Suppose that there are given limitations , Qf , 

Pi ! Pi ' C ifc' C ifc ■ 

(10) Qr<ft<Q, + , iel'; 

(11) Pi < Pi < pf, i^V; 

(12) c- k < Clk <c+, (i,k)eE; 

and the other restrictions can be represented by inequalities with given a~ k , af k 
for the functions a,ik{pi, Pk, <lik, Cik, dik) which have to be calculated: 



(13) a ik < a ik {pi,pk,qik,Cik,dik) < af k , (i,k) e E. 

The considering objective function depends both on flow and on potentials 
Pi,Pk, on intensities Qi, continuous Cik and discrete parameters dik - The problem 
consists in 

(14) minimize F = ^ F ik ' (Pi,Pk,qik,c ik ,d ik ) + F,- (2) (p t ,Qi) 

(Lk)eE i€V 



(15) subject to (El- EES]). 

The set of available values of discrete parameters dik in © means that a family 
of functions fd ik can act on the edge (i, k), and we have to select the best form 
of every equation ([5]). Thus, we have to select the best equation from the 
point of view of functional (|14p . We may interpret this as a selection of the 
most profitable equipment which is installed or can be installed on the place of 
the edge (i, k). 

The continuous parameters Cik in (|8]) can be interpreted as parameters that 
smoothly regulate the work of equipment dik in bounds (|12p . We may interpret 
(|13p as restrictions on the power, temperature, dissipation and other character- 
istics of the equipment that is represented on the edge (i, k) . The inequalities 
(|TU)l and the dependence of objective function on intensities Qi mean that the 
most profitable values of supplies and demands have to be chosen. 

The functions F^l , fd ik , a-ik staying in the objective function and in the re- 
strictions are non-linear and can be non-smooth and non-convex but expensive 
to compute. Usually the number of nodes and edges are hundreds with a ten- 
dency to be thousands. So the formulated problem is a problem of large-scale 
nonlinear discrete-continuous programming in general networks. This problem 
generalize the well-known minimal cost network flow problem. 
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3 Methods of the discrete-continuous nonlinear 
network optimization 

The formulated problem (0 - RH)) is a problem for searching a conditional cx- 
tremum for a function with continuous-discrete parameters of optimization. Ac- 
cording to the mixing character of variables, the combination of the continuous 
and discrete optimization methods is used for the problem solving. The integer 
and nonlinear programming and graph theory are basis of the offered method. 
Its main characteristic feature consists in the obtaining of dominant solutions 
on fragments of the network. 

3.1 Simple examples 

Let discrete parameters be fixed. The next 2 examples and orders are dual to 
each other. 

Example 1 Let a node v roo t is chosen. There exists a linear order £ p on V 
such that if we enumerate the nodes with respect to this linear order then: 1 ) 
the root has the number n Vroot = 1; 2) for every node with number n > 1 there 
is just one adjacent node with number n u < n. By other words, it is possible 
to construct an optimization process so that among potentials pi only p roo t is 
considered as an independent variable , and all other pi are dependent variables 
computing after p roo t and qik ■ 

Example 2 There exists a linear order £ q on nodes V such that if we enumerate 
the nodes with respect to this linear order then for every node with number n the 
numbers of all adjacent nodes except may be one are less than n : n u < n . By 
other words, it is possible to construct an optimization process so that among 
flows qik only non-tree flows are considered as independent variables, and all 
other flows qik are dependent variables. 

In terminology of linear programming, these independent variables corre- 
spond to non-basic variables, and .dependent variables are basic ones. 

3.2 The continuous dominant problems 

Let D and X are correspondingly the sets of discrete and continuous indepen- 
dent variables. Denote a problem for searching a feasible solution ([6] - fl3|) by 
P(D,X) and an optimization problem dU-E]) by FP(D,X). 

Problem P" dominates another one P and is called as a dominant for P, if 
the feasible set for contains such set for P. 

Lemma 1 For any discrete-continuous problem P(D, X) there exists a contin- 
uous dominant P"(Xo, X), where a set of discrete variables D is replaced by a 
set of continuous variables Xp. 

Proof. On the edges with a discrete variable dik, the equation tik — Pi — Pk 
replaces the family of equations (0 - ^ . A component tik of tension vector 
replaces the discrete variable dik and the continuous ones Cik there. H 

As consequence, there is the continuous dominant FP~ (Xd. X) for an initial 
optimization problem FP(D, X). 
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Assign a node where potential is either given or varied as a root. Let a set 
of discrete variables D is participated onto two classes D' and D c . Collections 
d! G D' are called as fragments of d e D. Consider an initial discrete-continuous 
problem ([5] -US]) by fixed d! € D'; denote it P d > (D c , X). We can construct for 
a problem Pd>{D c , X) its continuous dominant Pd'"(X D c,X). The problem 
Pd> " defines a feasibility of fragment d! . Fragment d! is infeasible, if P<j/ * has no 
solution. The explicit definition of fragment d! is based on the follows 

Lemma 2 Let G = (V, E) is a connected graph, vq € V is a node, and a finite 
set Mq C E consists of \Mq\ edges. Consider M as a set of subsets of Mq and 
write down d € M as 

(16) d= {di,...,d\ Mo \) p , dj=0,l. 

for a linear order p on Mq. There exists such a linear order it on Mq that for 
every m-length fragment d m = (di, d m , 0, 0)^ there is a connected subgraph 
Gi = (Vi,Ei),i = i(m), that is a minimal connected subgraph containing the 
node vq, all the edges from d m , and all their endnodes. 

Proof. A proof is constructive. It consists of the following steps: 

1) Construction a linear order £ v on a rooted spanning tree; 

2) Construction a linear order ir on the subset Mq of edges; 

3) Construction a minimal subgraph for a set d m of the first m edges from 
Mo - ■ 

Let the discrete variables are acting only on edges. Take the set of edges 
possessing discrete variables in the role of Mq. Let an edge e € Mq has a number 
j in a linear order tt; denote it as ej. With the notation (fT6"|) for d C Mq and 
d m = (d%, d m , 0, 0) for m - length fragment, let us allow to use a non-null 
value of a discrete variable acting on the edge ej G d as dj , and to use dj = for 
ej d . In the last case we shall say that discrete variable is not given there. So 
the linear order 7r on Mq induces a linear order on a set D of discrete variables, 
and we can tell about m-length fragment of discrete variables. The lemma [5] 
can be reformulated as follows. 

Theorem 3 Let G = (V,E,v roo t) is connected network with one node marked 
as v root . Let some edges of the network possess discrete variables, and the last 
ones are only there. There exists a linear order on the set of discrete variables 
D that: 

1 ) every d G D can be written as 

(17) d=(d!,...,du), di = 0,1,..., N h 

where di — means that value of corresponding variable is not given; 

2) for every m-length fragment d m = (di, d m , 0, 0) there exists a con- 
nected sub-network Gi — (Vi,Ei),i = i(m), that is a minimal sub-network con- 
taining the root, all the edges with non-null discrete variables from d m , and all 
endnodes of these edges. 

The solving method for the continuous-discrete problem FP(D, X) can be 
presented as a branching multi-level computational process. On the main branch 
of this process we solve the continuous problem FP~(D,X) or P{D,X). 
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3.3 A method for solving the continuous problems in gen- 
eral networks 

In continuous problems FP(X), FP~ (Xjj, X), and P<j'" j the components of 
continuous variable vector x — (xi, x n ) are divided into 3 groups: those from 
the initial problem FP(D, X); the flows in chords; and the continuous variables 
which change the discrete ones in the dominant continuous search. The flows 
in chords are of a special meaning. They can be computed as a solution of the 
equation system ^ , © , (JSJ by the fixed both discrete and continuous variables 
of 2 other groups. Then the dimension of continuous variables and of equation 
systems ([!]), (0) can be reduced thanks to the direct solving of such systems on 
a tree. 

To solve the continuous search problem, the classical methods of constrained 
and unconstrained optimization from nonlinear programming are implemented. 
The preference is given to the methods without derivatives because of non- 
convex and non-smooth functions in restrictions and objectives. 

Components of objective, penalty and Lagrangian functions are computed 
as soon as a network state component Si is been known, s = (qik, Qi,Pi, Cjfc, dik). 

By computation of state s, we have to distribute flow and potential. To 
effectively solve the equation systems ^ and (JSJ) on tree, respectively the width- 
and depth-first searches in graph are used. 



3.4 Modification of the branch-and-bound method 



In our modification of the branch-and-bound method, only m-length fragments 
are used. If fragment is infeasible then all its extension (d±, d m , j m +i, Jm) 
are excluded from the further examination. The examination is gone in a lex- 
icographic order of vectors d — (di, dj^)- Starting from d^ — (0,...,0), a 
transition from d^ to d^ n+1 ^ is performed by the next rule: if dy 1 ' = d m is 
fragment with length m < M and dy 1 ' is feasible then 



(18) 

else d( n+1 ) 
(19) 



= (d 1 ,...,d m) l,0,...,0) 
d m e m . Here e m = (0, 1, ...0) and 
d m ^i,d m + 1,0,..., 0), 



(di, 
(di, 



,dp_i + l,0 J ...,0), 



(0,...,0), 



if d m < N m ; 
if d m = N rn and 
3/i : 1 < fx < m : 
d^-i < Np-i & d k = N k , 
k = fi. . . m; 
tfdk=N k , Vfe=l,...,M. 



3.5 Another procedure for searching an optimal integer- 
feasible solution 

The idea of another procedure for FP(D, X) is connected with a possibility to 
check, if there is a realization of equation Uk = Pi — Pk in the equation family 
dS]) ;((n]) ■ This sub-problem is to construct mapping 
(20) 

{(tikiPi,Pk,Qik) I Uk =Pi ~Pk} — > {(dtk,c ik ,pi,p k ,,qik) | ©, ©, (H2])} . 
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The second procedure for solving FP(D,X) is proposed as follows. The dom- 
inant problem FP~ (Xjj, X) is to be solved. On the every step of continuous 
optimization FP" (Xd,X), there is checked the existence of feasible discrete 
variable dik and perhaps continuous ones Cik for every continuous variable Uk 
changing dik in the dominant problem. If there is no such feasible discrete vari- 
able dik then the corresponding continuous variable f,& would be penalized. In 
other words, Uk is penalized if there is no mapping (|20p . Then an optimum of 
problem FP~ (Xjj, X) brings an optimum for the initial problem FP(D, X). 

The difference between this method and the branch-and-bound approach is 
clear for the cases of pure discrete problems on a tree. Such class of problems 
means that graph G = (V, E) is a tree and there is no continuous variables at 
all: X = 0. Then the formulated version of branch-and-bound method can be 
very efficient. It brings a global optimum as a solution. 

For the same problem FP(D) on a tree, the second method based on the 
checking of integer- feasibility by the solving of dominant problem FP"(Xjj) 
can fall in a local optimum. To avoid it, a restart can be used. By restart, it is 
worth to move the root into another node, if the restrictions (fTTj) allow it. 

However, for the network with complicated topology the branch-and-bound 
method could be expensive if not prohibited. This is valid especially for the 
networks with a lot of cycles. Then the second procedure with penalizing non- 
realizable to, is available there. In practice, the method avoids the local optima. 

4 Risk, sensitivity and stability analysis for re- 
liability study 

There are 2 extreme approaches by analysis either a network state is stable. 
The first one is used by operating when a control vector is chosen (fixed) .and 
it is necessary to check either this control vector still lead to a feasible network 
state if initial or boundary conditions are changing. The second one is used by 
design when it is necessary to check either there is a value of a control vector 
for every possible operating point. The first approach we shall call the strong 
stability while the second one we shall call the weak stability. 

Let be given a value ttq of analyzing parameter tt, a. fixed control Uo, and 
So > 0. The scalar parameter tt is called feasible strong stable in the point 
7To by the fixed control uq if there is an interval (tto — S,tto + 5) , 5 > 0, such 
that the network state x(uo,tt) is feasible for every tt G (tto — S, tto + S). If 

5 > Sq then the network state xo — x{uo,tto) is called feasible strong stable 
with respect to the scalar parameter tt by the fixed control uq- 

To prove the stability of parameter tt in the neighborhood of ttq , we find 
such interval [7ri,7r2] containing tto that the network state x(uo,it) is feasible 
for every tt G [7ri,7r 2 ] , and for every small e > there are tt\ G [tti — £,7Ti] , 
TT2 G [tt2, 7T2 + s] that for 7r = tc%, 7T2 the network state x(uq, 7t) is not feasible. 

Let be given a value ttq of analyzing parameter 7r, a neighborhood Uo of a 
control uo and So > 0, rj > 0. Let F (u) — F (u,x (u, 7r)) be an objective function. 
The scalar parameter tt is called feasible weak stable in the point tto if there 
is an interval (tto — S, tto + S) , S > 0, such that for every tt € (ttq — S, tto + 5) 
there is a control u G Uo which provides that the network state x(u, tt) is feasible 
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and that 



(21) 



F (uq, x(uo, 7r) — min F± (u, x (u, w)) 

ueUo 



< rj 



where F\ (u, x) = F {u, x) + ip (u) and ip ( u ) 1S a cos t of control switch from uq 
to u. 

For a scalar parameter tt, it is possible to compute the above interval [tti, 7^2] 
showing a complete feasibility interval with respect to parameter n. For a set 
of analyzing parameters n = (tti, n n ) , a statistic test method could be used. 
If there are Si > that for the most points ir G II (71^0 — Si,iTio + 5i) (saying, 
for 95 %) there is a control u 6 Uq providing feasibility of the network state 
x(u, 7r) and condition (l?T1) then network state is called weak stable with respect 
to parameter set ir. The definition of strong stability with respect to parameter 
set 7r is evident. 

In complement to statistic approach, one simple result for 'definite' approach 
has to be mentioned. The capacity q% of an edge (i, k) can be found from 
equation © varying p t ,pk, d ik : 

!qtk I fd ik (puPk,qik) = 0, 
d lk e {i,...,AT ifc }, 
p n e\p-,p+], n = i,k 

In some cases we can precise potential intervals \p~ , pf] in (jlip and capacities 
therefore. 

Lemma 4 Let be given intervals \p~ , pf~\ ill)), functions fd ik in models (0) are 
monotonous relative both potentials Pi,Pk and qtk separately, and every model 
fdik i n (W * s chosen: dik — Const, (i,k) G E, i.e. fd ik — f ik - For a known 
flow qik , it can be found such intervals 

(23) [pTSpT"] c [ P -,pf] 

that the inequality system £8\) A1 1\) has a solution if and only if the system 
(0); flSP has a solution. 

In other words, for a node j G V every p° ^ [p™ ln ,p™ ax ] leads to a solution 
of equation system {JEl),Pj = P®} which violates the inequality system ©, (ITT1) . 



5 Conclusions 

The solution methods and their implementation are described for a class of 
problems of large-scale nonlinear discrete-continuous network optimization in 
general networks. Objective function depends both on flow and potential. Min- 
imal cost and maximal flow problems are generalized. The problem consists 
both in the optimal choice of dependencies on flow and potential from the fam- 
ilies given on some edges and in optimization of values of intensities, flow, and 
potential. 

The formulations, methods and algorithms are developed for the follow- 
ing problems: 1) the (pure) continuous nonlinear optimization in general net- 
work with a given functional dependence between flow and potentials; 2) the 
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(pure) discrete nonlinear optimization problem of selection the best functional 
dependence between flow and potentials from given families; 3) the (mixed) 
continuous-discrete nonlinear optimization in general network with given fam- 
ilies of dependences between flow and potentials. The developed optimization 
method represents a branching multi-level computational process. It is based 
on nonlinear and integer programming and on graph theory. Its main character- 
istic feature consists in obtaining of dominant solutions on network fragments. 
During optimization, the offered methods adapt themselves to the structure of 
every network and to the features of both objective function and variables. It 
allows to use such objectives as minimization of cost, of expenditure, and of 
set-point deviation, maximization of flow, of profit, and so on. 

On the base of the developed optimization methods, some approaches for 
sensitivity analysis, stability investigation and reliability study are proposed. 
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